#include <iostream>

#include "Crible.hpp"
#include "Paire.hpp"
#include "Factor.hpp"
#include "Erathostene.hpp"
#include "Polynome.hpp"

using namespace std;


Factor * Crible::sieving(mpz_class inLowBound, mpz_class inUpBound, int inNbelementRationnal, Erathostene & erathostene, Polynome & poly) const
{
   
   int nbCouple = 0;
   
   
   mpz_class smooth = 0;
   
   mpz_class pgcd = 0;
   mpz_class indiceJ = 0;
   mpz_class b = 1;
   mpz_class cptK;
   mpz_class indiceMax = 2 * inUpBound + 1;
   mpz_class valeur;  //Utilie dans le calcul du for
   
   Paire paireCourante;


   
   //Calcul du nombre total d'élément des bases + 1 
   int nbMaxCouple = _algebraicFactorBase->getPaire() 
                   + _rationnalFactorBase->getPaire()
                   + inNbelementRationnal
                   + 1;
   
  
   Factor * listeCouple = new Factor(500);
   Paire tmp;
   
   mpz_class * tab1 = new mpz_class [mpz_get_ui(indiceMax.get_mpz_t())];
   mpz_class * tab2 = new mpz_class [mpz_get_ui(indiceMax.get_mpz_t())];
   
   mpz_class * courant;
   mpz_class * courant2;
   
   mpz_class * fin2;
   mpz_class * fin;
   
   fin = &tab1[mpz_get_ui(indiceMax.get_mpz_t()) - 1];
   
   while(nbCouple < nbMaxCouple)
   {
      
      // set a[i] = i + bm for i ∈ [−C; C]
      valeur = inLowBound;
      
      for(courant = tab1 ; courant <= fin ; courant++)
      {
         *courant = valeur + b * _m;
         if(*courant < 0) *courant = (-1) * (*courant);
         valeur++;
      }
      
      for(int i = 0 ; i < _rationnalFactorBase->getTaille() ; i++)
      {
         //On recupere la paire courante
         // Sous la forme (p,r)  
         paireCourante = _rationnalFactorBase->getPaire(i);
         
         /*
         Divide out largest power of p from a[j] where
         j = −bm + kp for k ∈ Z so that −C ≤ j ≤ C.
         */
         mpz_class p = paireCourante.getPaire(0);
         
         smooth = erathostene.higherFactor(p);
//          cout << smooth << " " << p << endl;
         //Premiere boucle pour les valeurs negatives
         
         
         // Recherche du premier indice
         cptK = (inLowBound + b * _m);
         cptK /= p;
         
         
         // Calcul du premier indice
         indiceJ = (-1) * b * _m + cptK * p;
         if(indiceJ < inLowBound)
            indiceJ += p;
         
//          cout  << "1" << endl;
//          cout << indiceJ << endl;
         // Recherche de l'indice dans le tableau
         cptK = inLowBound;
         courant = tab1;
         while(cptK != indiceJ)
         {
            cptK++;
            courant++;
         }
         
         
         while(indiceJ <= inUpBound)
         {
            // Operation sur a
             if(*courant != 0)
             {
                while(mpz_divisible_p((*courant).get_mpz_t(), smooth.get_mpz_t()))

               {
                  *courant /= smooth;
               }
             }

            
            // Indice suivant
            indiceJ += p;
            while(cptK != indiceJ)
            {
               cptK++;
               courant++;
            }
         }

         
      }
      
      // set e[i] = bdeg(f ) f (i/b) for i ∈ [−C; C]
      valeur = inLowBound ;
      fin2 = &tab2[mpz_get_ui(indiceMax.get_mpz_t()) - 1];
      
      for(courant2 = tab2 ; courant2 <= fin2 ; courant2++)
      {
         *courant2 = poly.bpowdf(valeur, b);
         
         
         if(*courant2 < 0) *courant2 = (-1) * (*courant2);
         valeur += 1;
      }
      

      for(int i = 0 ; i < _algebraicFactorBase->getTaille() ; i++)
      {
         //On recupere la paire courante
         // Sous la forme (p,r) 
         /*Divide out largest power of p from e[j] where
         j = −br + kp for k ∈ Z so that −C ≤ j ≤ C. */
         paireCourante = _algebraicFactorBase->getPaire(i);
         
         mpz_class p = paireCourante.getPaire(0);
         mpz_class r = paireCourante.getPaire(1);
         
         smooth = erathostene.higherFactor(p);
         
         // Recherche du premier indice
         cptK = (inLowBound + b * r);
         cptK /= p;
         
         // Calcul du premier indice
         indiceJ = (-1) * b * r + cptK * p;
         if(indiceJ < inLowBound)
            indiceJ += p;
         
         // Recherche de l'indice dans le tableau
         cptK = inLowBound;
         courant2 = tab2;
         while(cptK != indiceJ)
         {
            cptK++;
            courant2++;
         }
         
         
         while(indiceJ <= inUpBound)
         {
            if(*courant2 != 0 )
            {
               while (mpz_divisible_p((*courant2).get_mpz_t(), smooth.get_mpz_t()))
               {
                     *courant2 /= smooth;
               }
            }

            
            // Indice suivant
            indiceJ += p;
            while(cptK != indiceJ)
            {
               cptK++;
               courant2++;
            }
         }

      }


      
      
      /* for i ∈ [−C; C] 
      i. if a[i] = e[i] = 1 and gcd(i, b) = 1 add (i, b) to rels  */
      
      valeur = inLowBound ;
      courant2 = tab2;
      
      for(courant = tab1 ; courant <= fin ; courant++)
      {
         mpz_gcd(pgcd.get_mpz_t(), valeur.get_mpz_t(), b.get_mpz_t());

//         cout << "(" << valeur << " ; " << b << ") a : " << *courant << " ; e " << *courant2 << endl;
         if( *courant == *courant2 && *courant == 1 && pgcd == 1)
         {
            tmp.setPaire(valeur, b);
            listeCouple->add(tmp);
            
             cout << "Couple : " << tmp << endl;

            
            ++nbCouple;
         }
         ++courant2;
         valeur += 1;
         
      }

      b++;
       cout << "nombre de couple " << nbCouple << " ; nombre max " << nbMaxCouple << endl;
   }
   
   
   
   delete [] tab1;
   delete [] tab2;
   
   return listeCouple;
   
}

